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Abstract. Motivated by the particle acceleration problem in pulsars, we numerically investigate electrostatic 
instability of electron-positron pairs injected in an external electric field. The electric field is expected to be so 
strong that we cannot neglect effects of spatial variation in the 0-th order distribution functions on the scale of 
the plasma oscillation. We assume that pairs are injected mono-energetically with 4- velocity uo > in a constant 
external electric field by which electrons (positrons) are accelerated (decelerated) . By solving linear perturbations 
of the field and distribution functions of pairs, we find a new type of electrostatic instability. The properties of the 
instability are characterized by uo and the ratio R of the braking time-scale (determined by the external electric 
field) to the time-scale of the plasma oscillation. The growth rate is as large as a few times the plasma frequency. 
We discuss the possibility that the excited waves prevent positrons from returning to the stellar surface. 
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1. Introducton 

Waves in homogeneous plasmas are well described by 
the linear perturbations that have the Fourier-harmonic 
dependence of the form exp [i(k ■ x — u)t)]. The properties 
of various wave modes have been extensively studied for 
various physical situations. Plasma instabilities, such as 
the tw o-stream instability, the Weibel instability l|Weibel L 
and many others have been recognized as impor- 
tant processes in astrophysics as well as in laboratory situ- 
ations. For inhomogeneous plasmas, the Fourier-harmonic 
dependence is not assured in a strict sense. If the wave- 
length is short enough compared with the scale of the 
inhomogeneities, we may neglect effects of the spatial gra- 
dient of the distribution functions of particles, and carry 
out the Fourier-harmonic expansion. We call this treat- 
ment the local approximation hereafter. 

When static electric field exists in a plasma, it accel- 
erates particles and leads to inhomogeneous velocity dis- 
tributions of particles. Wave properties in an electric field 
were studied in geophysical researches, adopting the lo- 
cal approximation (e.g. Misra and Singh 1977, Misra at 
al. 1979, Das and Singh 1982). In many high-energy as- 
tronomical phenomena, electric field is expected to be so 
strong that the local approximation is not adequately ap- 
plied. Such a situation typically appears in the magneto- 
sphere of pulsars. 

A spinning magnetized neutron star provides huge 
electric potential differences between different parts 



of its surface as a resu lt of unipolar induction 
l(Goldreich and Julian L Il969l) . A part of the potential dif- 
ference will be expended as an electric field along the mag- 
netic field somewhere in the magnetosphere. Although a 
fully consistent model for the pulsar magnetosphere has 
yet to be constructed, several promising models have been 
considered. Among them, the polar cap model l|Sturrock I . 
Il97li iRuderman and Sutherland"! Il97fj) assumes that an 
electric field E\\ parallel to the magnetic field lines exists 
just above the magnetic poles. The electric field acceler- 
ates charged particles up to TeV energies, and resultant 
curvature radiation from these particles produces copious 
electron-positron pairs through magnetic pair production. 
These pairs may provide gamma-ray emission by curva- 
ture radiation or synchrotron radiation as well as coherent 
radio emission and a source for the pulsar wind. 

The localized potential drop is maintained by a pair 
of anode and cathode regions. In the cathode region 
the space charge density p deviates from the Goldreich- 
Julian (GJ) density /?gj — —£IB z /2ttc negatively, where 
Q = 2tt/T is the angular velocity of the star and B z 
is the magnetic field strength along the rotation axis. 
On the other hand, p deviates positively for the anode. 
Outside the accelerator the electric field is screened out. 
In the polar ca p model, especially for space charge lim- 
ited flow model ifFawlev et al. Lll97'7llScharlemann et al. L 
ll978tlArons and Scharlemann Lll979h . where electrons can 
freely escape from the stellar surface, i.e., E\\ = on the 
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stellar surface, the formation mechanism of a static pair 
of anode and cathode, which can sustain enough poten- 
tial drop for pair production, is a long-standing issue. 
Current flows steadily along the magnetic field line so 
that the charge density is determined by the magnitude 
of the current and field geometry with suitable bound- 
ary conditi ons. Good exam ples for space charge limited 
flow are in Shibata (1997). When pqj < and the elec- 
tron density (n oc B/v ~ B/c, where B is the mag- 
netic field strength) is larger than the GJ number density 
(tjqj = |pGj/e| oc B z , where — e is the electronic charge) 
on the stellar surface, a cathode is provided on the stel- 
lar surface. The cathode accelerates electrons. When the 
field lines curve away from the rotation axis, n deviates 
riQj nagatively so much more for 'away' curvature, which 
enhances the cathode. Hence electrons continue to be ac- 
celerated, and potential drop becomes large enough to 
produce pairs. The mechanism of the electric field screen- 
ing, i.e., a way to provide an anode, has been considered 
to be provided by pair polarization. Although most pa- 
pers take it for granted that copio us pair production can 
insta ntly screen the field, recently IShibata et al~l fl998, 
2002) castcd doubt on this issue; the electric field screen- 
ing is_jiat_an_easy task as cons idered usually. 

IShibata et al. I (^398, 200^) investigated the screening 
of electric fields in the pair production region. They found 
that the thickness of the screening layer is restricted to be 
as small as the braking distance Ze = m c c 2 /\eE\\ | for which 
decelerating particles become non-relativistic, where m e is 
the electron mass. If the above condition does not hold, too 
many positrons are reflected back and destroy the nega- 
tive charge region (cathode). In order to screen the electric 
field consistently, huge number of pairs should be injected 
within the small thickness Ze- The required pair multipli- 
cation factor per one primary particle is enormously large 
and cannot be realized in the conventional pair creation 
models. Thus, some other ingredients are required for the 
electric field screening. 

In the previous studies of the screening, pairs were 
assumed to accelerate or decelerate along the 0-th or- 
der trajectories determined by Eu . However, if an electro- 
static (longitudinal) instability occurs, the excited waves 
may produce effective friction. Friction on particles change 
the charge polarization process. Thus, instability in the 
presence of an external electric field may have a rel- 
evance to this problem, which motivates us to make 
an exploratory study in this paper. Various instability 
mechanisms outside the accelerator have been studied 
for pair plasmas along with a primary b eam in r e lation 



to coherent radio emission mechansims llHinata 
Cheng and Ruderman L Il977t lAsseo et al. 



1976; 
1983; 



is Ze — 10 2 cm, while the length-s cale of plasma oscilla- 
tion is (TO c c 2 /47re 2 n G j) 1/2 -lcm ijShibata et al. L IT9981 
120021) . where uqj — \pQj/e\. Particles are accelerated or 
decelerated in a period that is shorter than the typical 
time scale of plasma oscillation. Therefore, the distribu- 
tion function is not uniform on the scale we consider. The 
local approximation is not adequate to deal with plasma 
oscillation. 

Properties of a pair plasma in such a strong electric 
field have not been studied. Such studies may bring us 
a new key to understanding astrophysical phenomena. 
One of our purposes is to examine if electrostatic insta- 
bility makes the screening easier. As the first step to- 
ward this purpose, in this paper we investigate electro- 
static instability of pairs injected in an external electric 
field. Investigations of electrostatic waves, when we can- 
not adopt the local approximation, may be important not 
only for pulsars but also for other high-energy astronom- 
ical phenomena. Since an analytical treatment is difficult 
in this case, we simulate electrostatic waves numerically 
in idealized situations. In §2 we mention the two-stream 
instability without an electric field for comparison. In §3 
we describe the situation we consider. The most simpli- 
fied physical condition is adopted; the electric field, injec- 
tion rate, injection energy are constant. In §4 we explain 
our numerical method. Our method can treat only linear 
waves. Numerical results are summarized in §5. Our re- 
sults show a new type of plasma instability due to electric 
field. §6 is devoted to summary and discussion. 

2. Two stream instability with local approximation 

First of all, for reference, we consider one-dimensional (1- 
D) homogeneous flows of electrons and positrons in the 
absence of electric field. The pair-distributions are func- 
tions of the 4- velocity u — 0/ 'yT— /3 2 , where /3 = v/c. 
For simplicity, we assume that the distribution functions 
are expressed by the step function ©(it) as 

/_ = C [0(it — it 3 ) — 0(u — 114)] for electrons, (1) 
/+ = C [Q(u — ui) — <d(u — U2)] for positrons, (2) 

where 114 > U3 > 112 > u\ with constant U{. Electrons 
homogeneously distribute between U3 and 1/4 in it-space, 
while positrons similarly distribute between iti and 112 (see 
Figure 1). As will be seen in the next section, the distribu- 
tions in an external electric field that we will consider are 
similar to the above distribution functions. In the linear 
perturbation theory, the dispersio n relation for electro- 
static waves l)Baldwin et al. L Il9691 is given by 



Lvubarskii U l992'l. However, plasma instability inside the D(k,Lo) 1 
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accelerator has not been studied. The Lorentz factor of the 
primary beam (r ~ 10 6 -10 7 ) is much larger than that of 
electron-positron pairs (7 ~ 10 2 -10 3 ). In such a case, it is 
difficult to induce the two-stream instability. However, it is 
not clear whether pairs stably flow in the electric field E\\ 
or not. For typical pulsar parameters, the braking distance 



(ft 



(3 df a 
u — kv du 



= 0, 



(3) 



where q a and m a denote charge and mass of the particle 
species a. Solutions of the dispersion relation usually yield 
a complex frequency. The imaginary part of ui corresponds 
to the growth rate of waves, oji. A positive growth rate 
Ui > implies an exponentially growing wave, while a 
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negative u>i does an exponentially damped wave. Adopting 
the distribution functions (TjJ an d ©, we obtain 

Q + " Wi) ~ /VP - kfo) + /3 3 /(w - fcft) 

-(3 4 /(u> - kfa) = 0, (4) 

where the subscripts i denote corresponding values to Uj, 
v = uj/ujp, k — kc/ujp, and uj 2 = Ane 2 C '/m c . The above 
equation is reduced to a quartic one. Although it is messy 
to obtain the solutions analytically, we can numerically 
confirm that uj has complex values, for a given k. The solu- 
tions obey the usual two-stream instability properties; the 
larger the difference in velocities of electron and positron 
flows is, the smaller the maximum growth rate and the 
corresponding k become. 

When u 2 =1*3, electrons and positrons distribute con- 
tinuously, though the average velocity is different. In this 
case, the dispersion relation is reduced to a quadratic 
equation, and we obtain 



Q = k{fii + P±)/2 ± - sJk?{p A - /3i) 2 + 4(/3 4 - fa). (5) 

Apparently, u) is real. Therefore, even though the veloc- 
ities of the two flows are different, instabilities are not 
excited in homogeneous pair plasmas, as long as electrons 
and positrons distribute continuously in it-space. This is 
because the absolute values of their charges and masses 
are the same for electron-positron pairs. As will be shown 
in the next section, the distribution of pairs injected in 
an external electric field is similar to the above distribu- 
tion locally. However, we will show that an electric field 
induces instability. 



3. Pair injection in an external electric field 

In a strong magnetic field as in the pulsar polar caps, 
transverse momenta of relativistic electrons and positrons 
are lost during a very short time via synchrotron radia- 
tion. These particles move along the magnetic field lines 
and their distribution functions are spatially I-D. In this 
section we consider I-D distribution functions of electron- 
positron pairs injected in an external electric field that is 
parallel to magnetic field lines. As the conventional the- 
ories for the two-stream instability implicitly assume, we 
neglect the toroidal magnetic field due to the global cur- 
rent of plasma. Only within this treatment the 1-D ap- 
proximation is adequate. In order to simplify the situation, 
we assume the external electric field Eq is constant. In pul- 
sar models, there exists a primary beam which produces 
current flow. The external electric field is determined by 
the complicated combination of the beam current, injected 
pair plasma, and GJ density. The constant E requires 
that the charge density of the beam, pairs, and GJ den- 
sity cancels out, which may be an artificial situation. As 
will be discussed in §6, the approximation of constant Eq 
is justified for a smaller rate of pair injection, which may 
be realized for actual pulsar parameters. 



Anyway, we depart from actual pulsar physics, and 
deal with plasma physics in an idealized situation here- 
after. We adopt the 1-D approximation and assume the 
existence of the background charge which leads to con- 
stant Eo. In our treatment we totally neglects effects of 
the existence of the background on development of waves, 
and consider the behaviour of pair plasma only. Pairs are 
assumed to be injected between z = and zm at a constant 
rate no- In our calculation the pair injection is monoen- 
ergetic with 4-velocity u — uq > 0. Let us start from as- 
suming the steady state of flows of electron-positron pairs. 
The distribution function fo(z, u) satisfies the Boltzmann- 
Vlasov equation 



df gE df 
dz rripC du 



h S{u - u Q ), 



(6) 



for < z < zm- Hereafter we assume Eq < for definite- 
ness. 

Then, injected positrons (q — e > 0) will be deceler- 
ated as 



Uq - (t - iinj)/*E 



(7) 



where tg = m e c/\eEo\, and i; n j is the injection time. The 
value of £e represents the time-scale in which the Lorentz 
factor of decelerating particles decreases by one. Positrons 
will be turned back following the trajectory 



z/Ie = 7o - 7 + Zinj/fe, 



(8) 



where 7 = yjv? + 1, 70 = ^/u^ + l, l E = c*e, and z in j is 
the injection point. 

Electrons (q = — e < 0) will continue to be accelerated 

as 



U = U +(t- iinj)/*E, 
Z/Ie = 7 - 70 + Zinj/lE- 



(9) 
(10) 



In our model electrons and positrons distribute in the 
2-D phase space (spatially 1-D) as is illustrated in Figure 
2. In the phase space the trajectory of electrons injected at 
z = provides the maximum 4-velocity % = «m(z) at z, 
while u of positrons injected at z — zm corresponds to the 
minimum 4-velocity u m = u m (z) at z. The region enclosed 
by the curves in Figure 2 is composed from the family of 
trajectories of pairs. As long as pairs are uniformly in- 
jected at a constant rate, the distribution function is con- 
stant because of the Liouville's theorem. The distribution 
functions between z = z rct = (70 — 1)^e (see Figure 2) and 
zm are the same as in Figure 1 with u\ = u m , 114 — mm, 
and U2 — U3 — uq. In this region the distribution functions 
are expressed as 

/o- = no [Q(u — uq) — Q(u — um(z))] for electrons(ll) 
/o+ = n [Q(u — u m (z j) — Q(u — uq)] for positrons(12) 

where hq = t-^hQ. It is easily confirmed that /q- and 
/o+ satisfy the Boltzmann-Vlasov equation JBJ). Electrons 
and positrons distribute continuously. The charge density 
given by Eqs. (11) and (12) is not strictly zero so that 
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our approximation of constant electric field is not self- 
consistent, unless the background charge cancels the total 
charge. If the background does not play such a role, our 
treatment is correct only when uq is small enough. The 
quantitative condition for n will be discussed in §6. As 
was discussed in §2, if we neglect the electric field and 
adopt the local approximation, no wave instability is gen- 
erated for this continuous distribution . However, in our 
consideration the time-scale ts is too short (or Ie is too 
short) to adopt the local approximation. 

In the region < z < z rot the pair distribution has 
separate two streams. Since this region is peculiar in our 
idealized model, we do not consider the waves in this re- 
gion hereafter. 

4. Numerical method 

We consider linear perturbations of the distribution func- 
tion and electric field as 



/ = h(z,u) + fx(z,u,t), 
E z = Eo+Ehfat), 



(13) 
(14) 



where |/i| <§; f and \Ei\ <C \E \. Since the unperturbed 
distributions of pairs are inhomogeneous, we cannot carry 
out a Fourier-harmonic expansion of the perturbations. 
Therefore, we solve time development of the perturbations 
rather than obtain the linear modes. As we have seen in 
the previous section, our idealized situations lead to the 
simple distribution function fo, which makes numerical 
computaions easier. 

We directly solve the perturbations fi and E\ from 
the linearly perturbed Boltzmann-Vlasov equation 



Dfx 
Dt 



dfi 
dt 



dfi , qE dfi 



dz 



m e c du 



qEi df Q 
m c c du 



(15) 



where Dt is the differential along the 0-th order trajectory. 
The Ampere-Maxwell law is written by 



dEi 



-47TJ, 



(16) 



where j is the perturbed current density. Initially E\ de- 
pends on only z, and the magnetic field B x = B y = 0. 
Then, the Faraday law ensures that B x and B y remain 
zero all the time. The component B z does not affect time 
developments of E z and f±. 

We set up grids along the 0-th order trajectories of 
pairs in the 2-D phase space (s = z/l-£,u). Following the 
Lagrangian method we follow time evolution of fx in these 
grids from equation (|15fl . For s rct < s < sm, equation (|15|l 
is rewritten as 



DF- 

Dt 
DF+ 

Dt 



= E(s) [S(u - u ) - 8(u - uu{s))\ , 
= -£{s) [8{u - u m (s)) - 5(u - m )] 



(17) 
(18) 



where we introduce dimensionless values F = /i/e, t = 
t/t-E, and £ = l-E,noEi/\Eo\. At the injection of pairs 
(it = uq) F obtains a finite value in proportion to the 



perturbed electric field E. Then, the value of F is propa- 
gated along the 0-th order trajectory (see Figure 2), i.e., 
F is conserved along the characteristics. The disturbance 
F + initially propagates forward, and then turns back at a 
distance 70 — 1 from the injection point, while F_ simply 
propagates forward. For \u\ < uq, F + (u,s) is originated 
from inner injection positions between s — 70 + 1 and s. 
On the other hand, for \u\ > Uo, F+(u,s) is originated 
from the outer injection positions > s. The values of F on 
the trajectories of the pairs injected at z = and z = zm 
are special and they will continue to be changed by the 
electric field during propagation. 

On the other hand, the evolution of electric field is 
calculated by the Eulerian method; 



9E „ 2 



duf3-F_ 



du(J+F+ 



(19) 



where R = uj p t-E with the plasma frequency lu p = 
y / 47re 2 no/m e . The parameter R is the ratio of the braking 
time-scale to the typical time-scale of the plasma oscilla- 
tion. In our case R cx t^ 2 oc |i?o| -3 ^ 2 - In the pulsar polar 
cap model electric field is so strong that R is much smaller 
than one. Since the charge conservation is assured by the 
Boltzmann-Vlasov equation, the Gauss law, 



as 
a7 



R' 



duF, 



duF_ 



(20) 



is automatically satisfied, if the law is initially satisfied. 
The configuration of E is determined by the charge density 
that is an integral of F + —F_ . Since the propagations of F + 
and F_ are complex, it is difficult to predict the behaviour 
of £ intuitively. 

We have ascertained that results obtained from our 
numerical code satisfy the Gauss law. In addition we have 
checked our code by reproducing two stream instability in 
the absence of electric field, using the distribution func- 
tions in §2. 

5. Results 

We have simulated electrostatic waves from various ini- 
tial conditions and parameter values. We are interested 
in a parameter region R < 1. In this region the typical 
wavelength of plasma oscillation ~ l p = c/uj p is longer 
than the braking distance Ie = Rip- We give an initial 
disturbance in a spatially limited region. As will be shown 
below, when Ruq >^ 1, we find an absolute instability 
in which disturbance grows in amplitude but always em- 
braces the original region, where the initial disturbances 
of F and £ are given. The condition Ruq > 1 means that 
the distance injected positrons move forward before they 
turn back, Ze(7o — 1) ~ 'e u 0i is larger than the length- 
scale of the plasma oscillation / p . On the other hand, for 
~ 0.1 < Ruq <~ 1, we find a convective instability in 
which disturbance grows while propagating away from the 
original region. The waves excited from the convective in- 
stability propagate backward. Empirically, the results do 
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not largely depend on the spatial size of the pair injection 
region Sm = zm/Ie, which determines the minimum 4- 
velocity u m (z). In this section we show some examples of 
the instabilities found in our simulations. The parameters 
and initial conditions are summarized in Table 1. 

The initial conditions are taken to satisfy the Gauss 
law. Given the parameters R, uq, and Sm, we set the initial 
values of the disturbances F and £ as 

£(s) = £(cos(fci(*-*i))-l), (21) 
F_(s,u) = {ki ~^ )£i S m(h(s- Si ))6(u-uo), (22) 

*+(*>«) - sin (h(s- Si ))5(u-u ), (23) 

for Sj < s < S[ + 2n/ki. 

The initial disturbances are confined within a small 
spatial region of one wavelength (27r/fc; <C In the 

other region there is no disturbance. The initial pertur- 
bation of £(s) has a single sign with a form of a cosine 
curve. On the other hand, F- and F + have the form of 
a sine curve, satisfying the Gauss law. The parameter &a 
induces asymmetry of the charge density of electrons and 
positrons. The total charge density (oc F + — F_) docs not 
depend on &a , and also has a form of a sine curve with ki . 
We have tried various values of k& in our simulations and 
find that the ratio of E_ to F + is settled as time passes 
irrespective of /ca- When instabilities occur, growing wave 
modes end up dominating other modes. Therefore, results 
do not largely depend on the initial conditions. 

First we describe the results for R = 0.1 (the calcu- 
lations RUN1-RUN3) and see the behaviour of the linear 
perturbations for various values of the parameter uq. In 
Figures 3 and 4 we plot electric field £ for RUN1 for which 
Ruo = 1. In this calculation, the initial disturbance exists 
from s; = 120 to s; + 27r/fc; ~ 180. Since positrons will 
turn the direction of motion after r ~ 10, we must follow 
the disturbance much longer than that. As is illustrated 
in Figure 3, at r — 25 the disturbance £ remains in the 
originally disturbed region. As time passes, the amplitude 
around the original region of the disturbance grows, and 
the wave packet spreads backward little by little. In the 
forward region s >^ 200, we do not observe any grow- 
ing wave. Although particles move almost at the light 
velocity, the disturbances remain around the original re- 
gion and the wave packet docs not spread at the light 
velocity. In order to show the growth of the amplitude, 
in Figure 5 we plot the time evolution of Em that is the 
electric field for the maximum amplitude. The maximum 
electic field Em oscillates over positive and negative re- 
gions. Initially Em changes complicatedly because of the 
initial conditions we artificially set. As time passes, -Em 
smoothly grows while oscillating. The period of the oscil- 
lation of Em is — 20£e = 2/w p . The growing time U, where 
| .Em | ~ exp(i/ij), is ~ 50i E = 5/w p . 

Let us look into the the behaviour on a shorter time 
scale for RUN1. At a fixed position s, the local elec- 
tric field £(s) grows while oscillating. However, when we 
see spatio-temporal behaviour, we notice that the spatial 



pattern propagates backward while changing their am- 
plitude. As is shown in Figure 4, waves propagate from 
s ~ 200 ~ Si + In fk\ backward with a growing amplitude. 
The amplitude becomes maximum around s; = 120, and 
then the amplitude declines while propagating backward. 
This decline leads to the confinement of the wave packet. 
Even though waves pass the disturbed region many times, 
waves exist only in a spatially limited region. In the wave 
packet of £(s) there are multiple peaks and bottoms, and 
the most prominent one of them corresponds to -Em- If 
we define the 'phase velocity' as the velocity of peaks (or 
bottoms), the phase velocity (~ —2.8c) turns out to be 
faster than the velocity of light. As peaks propagate back- 
ward, the peak or bottom associated with -Em alternates 
one after another, so that the position of Em hangs around 
the original region. When we define the group velocity by 
averaging the velocity of the position of £ = Em for a 
longer time scale than the oscillation period, the group 
velocity turns out to be almost zero. The wavelength A is 
about 60^e — 2-7rZ p which is almost the same as the initial 
wavelength of the disturbance. Even if we start from an- 
other ki, the growing mode dominates others and the final 
wavelength is the same as this result, 60Ze- The final wave- 
length of growing waves is unchanged for different initial 
conditions. Figure 6 shows charge density distributions at 
t = 195. The number density of electrons (positrons) has 
opposite (same) sign of the charge density in Figure 6. The 
phases of number densities of electrons and positrons are 
the same. The amplitude of the positron density is always 
larger than that of the electron density. The difference in 
the number densities of positrons and electrons is propor- 
tional to the total charge density which satisfies the Gauss 
law. 

Next we discuss the results of RUN2 (Ruo = 30 > 1). 
The initial distubancc ranges from s; = 500 to s; + 27r/A:i ~ 
560. In Figure 7 we show electric fields at several epochs. 
The wave profiles are not so simple compared to the case 
of RUN1 (Ruq = 1). There are waves propagating both 
forward and backward. This may be because the distance 
positrons move forward before they returns (~ 300Ze from 
their injection point) is longer than the typical wave- 
length 2irl p ~ 60Ze- Though the properties of the waves 
are complex, we can see that there is an instability in 
this case, too. The waves around the original region grow 
while diffusing both forward and backward. We note that 
a separate component of the disturbance appears around 
z < z rot (s < 300 in this case). This disturbance may be 
due to two stream instability, which grows faster than in 
the other region. 

Although we do not show here, for a much larger value 
of Ruo, absolute instabilities are found to occur in our 
simulation. However, we do not further pursue this issue 
because we need to calculate over a much wider region of 
s and resultant memory in computation becomes large for 
a large value of uq. 

In cases for Ruo = 0.3 < 1 (the calculations RUN3), 
absolute instability does not occur, but convective in- 
stability propagating backward occurs (see Figure 8). In 
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RUN3 the initial disturbance is given from s = 420 to 
~ 426 with a wavelength 27rlE/h ~ 6/e- The disturbance 
of the electric field is seen to propagate backward, and 
before long characteristic waves grow. The wavelength 
of the growing wave (~ 9/e) is slightly longer than the 
initial length and almost constant. As the disturbance 
propagates, the wave packet spreads and the number of 
waves in the packet increases. The growing time of Em is 
about IOO^e = 10/cjp (see Figure 9). The phase velocity, 
w p h — —1.09c, is faster than the velocity of light. The am- 
plitude of each peak in the packet initially grows. As the 
peak approaches the head of the wave packet, the growth 
of the amplitude turns over to damping. This behaviour 
is similar to that in the absolute instability for RUN1. We 
obtain the group velocity v g ~ —0.85c. Figure 10 shows 
that the charge density is dominated by positrons. We 
have simulated for Ruo = 0.3 with various initial condi- 
tions. However, in any case there is no sign of wave insta- 
bility propagating forward. 

We have investigated for R = 0.01 also (the calcula- 
tions RUN4-RUN6) , and confirmed that the qualitative 
results are determined by the value of Ru n (see Table 
2). The absolute instability and convective instability are 
induced for Ruq = 1 and Ruo = 0.3, respectively. For 
Ruq = 0.03 (RUN6), however, we do not find any insta- 
bility. A smaller value of R means lower density of pair 
plasma for a given value of the external electric field. In 
such a low density plasma, particles tend to be less affected 
by forces from other particles compared to the external 
electric field. Therefore, too small value of R <C 1/ito 
makes the plasma stable. From the above results, we 
may conclude that electrostatic instability is mainly de- 
termined by the parameter Ruq. The wavelength and the 
growing time are a few or ten times the typical scales of 
the plasma oscillation, Z p and l/w P , respectively. 

The value R has been assumed to be smaller than 1 
heretofore. We have also simulated for R > 1 and found 
instabilities. The wave properties are as complicated as the 
example in RUN2, so that we do not report the details of 
the results in this paper. When R > 1, Eq is so small that 

> l P - In the limit of E — (R — > oo), the instability 
does not occur as was shown in §2. However, we have not 
tried the cases R ^> 1, because of poor computational 
capacity so far. 

Judging from the backward spread of wave packets 
and the negative phase velocity, it is seen that return- 
ing positrons play a decisive role in the instabilities. The 
dominance of positrons in the charge density in compari- 
son with electrons also suggests that the instabilities are 
due to returning positrons. It is remarkable that positrons 
pass the same region twice, forward and backward. The 
typical wavelength of plasma oscillation (2ttI p ) may res- 
onate with the distance positrons move forward, (70 — 1)/e- 
As we have mentioned in §4, the excited electric field £ 
generates the disturbances of pairs, F(uo), at their injec- 
tion. The value of F(u) is transported along the trajec- 
tories of pairs, conserving their value. Since F + and F_ 
acquire the same value at their injection, the contribu- 



tion to the charge density is almost canceled out as long 
as electrons and positrons move forward together at ~ c. 
When positrons turn around, the charge is polarized. As 
for positrons, F + (u,s) at r conserves the information on 
the electric field £{s') in the past (r' = r — (u — u)), 
where s = s' + 70 — 7. The displacement between s and 
s' conforms to the trajectories of positrons. Therefore, the 
charge density of positrons (except for u = u m ) is pro- 
portional to a superposition of displaced £{s') at different 
times. This superposition will increase the amplitude of 
the charge density in response to evolution of £ . 

Since F + (u ^ u m ) is constant along the characteristics, 
the value F+ remains to be finite even after positrons enter 
the region where £ ~ 0, while F + (u m ) changes as long as 
£ exists. Our numerical results imply that F + (u m ) cancels 
out the charge density due to F + (u ^ u m ) in the region 
where £ ~ 0. This process prevents the waves from spread- 
ing backward at the speed of light. The propagation of the 
disturbance by electrons is relatively simple. The contri- 
bution due to F_(um) prevents the waves from spreading 
forward, as F + (u m ) does. In RUN3 and RUN5 (Ru < 1), 
positrons turn back and move backward quickly before 
F+(u m ) cancels out the charge density, so that growing 
waves propagate backward. 

6. Conclusions and Discussion 

In this paper we have found a new type of instability in 
electron-positron pair flows injected in an external elec- 
tric field, which is assumed to be spatially constant. The 
properties of the instability are characterized by the ra- 
tio R (the braking time-scale to the typical oscillation 
time-scale of the plasma l/w p ) and 4- velocity uo at injec- 
tion. For Ruq >~ 1 absolute instability is induced, while 
convective instability propagating backward is excited for 
^0.1 < Ruq <^ 1. The growing time in amplitude is 
as short as a few times the time-scale l/w p . The wave- 
length is also several times l p — c/uj p . The instabilities 
are caused by returning positrons. For Ruq <C 1, the pair 
plasma turns out to be stable. A small value of R implies 
that the plasma density is so low compared to the electric 
field Eo. For R <C 1/uo the collective interaction of the 
pair plasma is not important, so that each particle moves 
along the trajectory determined by E independently of 
other particles. 

Growing electrostatic waves may work as frictional 
forces. In this paper we have treated waves as linear per- 
turbations, following the propagation of disturbances in 
the distribution function and electric field. Our method 
does not allow us to follow processes of gaining or losing 
kinetic energy of each particle from the waves. The quasi- 
linear theory is not applied to deal with the reaction of 
particles as it is, because the disturbances do not have 
the Fourier-harmonic dependence. Thus, we consider the 
qualitative character of the effective reaction force from a 
numerical treatment as follows. The spatial averages of F 
and £ oscillate with time in our simulations. Therefore, 
the expectation values of F and £ can be considered to be 
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zero. On the other hand, when waves grow or are attenu- 
ated, the spatial average of the cross term, (F£ ) , may have 
a finite value. As is the case with the quasi-linear theory, 
the 0-th distribution function may change, following the 
2-nd order order approximations of the Boltzmann-Vlasov 
equation: 

d(f ) = fExdh 8G(u)^ 
dt m e c du du ' 

where 

G{u) = y-\ I dsF(u,s)E(s), (25) 
I y I J 

and (/o) = Uq independently of u in our case. In view of 
the Fokker-Planck approximation, G(u) oc {u)fo, where 
(ii) is the average change of u due to the reaction force 
per unit time. Since fo is constant in our simulation, G(u) 
is proportional to the reaction force. 

We plot G{u) in Figure 11 for r = 195 in RUN1. The 
modulation pattern of G(u) does not change, but the am- 
plitude grows with time. Apparently, the modulation pat- 
tern is asymmetric for particles of u > uq — 10 (electrons) 
and u < uq (positrons). These qualitative behaviour is 
common for the other RUNs. As Figure 2 and equations 
Ijl7(l and Ijl8(l show, perturbations are generated at u = uq 
and u = um (or u m ). Figure 11 shows G(u) for a region 
around u — uq only, and outside of this region G(u) has 
also significant value due to the disturbances generated at 
u = um and u = u m . However, the modulation pattern of 
G(u) for such regions oscillates with time. 

In the usual two-stream instability, the excited waves 
accelerate background fluid, and decelerate beam fluid. 
As is shown in Figure 11, the direction of the reaction 
force depends on u even in the same species of particles. 
The amplitude of G(u) takes always the maximum value 
at u = uq. If the effective reaction force grows enough, 
positrons (electrons) just injected (u ~ uq) feel positive 
(negative) force as is shown in Figure 11. Thus, the re- 
action force may make particles tend to stay around the 
regions of u = uq- The integral of G(u) around u — uq 
(roughly from u = —100 to 10 for Figure 11) for positrons 
is also positive. Such positrons are accelerated by the reac- 
tion force on average. However, the integral becomes neg- 
ative all the time, if we include the contribution due to the 
disturbances generated at u — u m , though G(u) for a large 
| it | oscillates with times. The returning positrons injected 
at z = zm feel a negative reaction force on average. On the 
other hand, the absolute value of the integral for electrons 
is much smaller than those for positrons. Therefore, the 
reaction force does not work as usual 'frictionaP force be- 
tween electrons and positrons. The particles just injected 
are most affected by the reaction force, and lose (or gain) 
their energy owing to waves. If the reaction force is strong 
enough, positrons just injected may have difficulty to turn 
back. If positrons suffer from such factional force, the dis- 
tribution function fo should be largely altered. This may 
help to solve the problem of the electric field screening in 
the pulsar polar caps. 



If excited waves grow enough to change trajectories of 
particles, the waves cannot be treated as the linear per- 
turbations. In a strong electric field, |/i| can be as large 
as fo before |£a| ~ \E \ achieves. We may simplify the 
energy-loss process of particles due to perturbed electric 
field as follows; particles from u = uq to uo + Au lose their 
energy owing to the waves at the same rate, where Au 
is the equivalent width of particles interacting the waves. 
Then, the growth rate of the field energy is roughly consid- 
ered as the average energy loss rate of pairs, which means 
tUiEf/4ir ~ |t Au\nom c c 2 , where 7 is the temporal change 
of 7 due to the reaction force. Here, we assume that the 
energy density of the waves attributed to induced par- 
ticle motions is comparable or negligible to the energy 
density of the electric field E\. When |7| > 1/£e, the fric- 
tional force is sufficient to alter trajectories of pairs. The 
above condition is rewritten as uiiEf/tu 2 > m c c\Eo\Au/e. 
Assuming ~ uj p ~ 10 8 , Au = 50, and |jEo| — 10 5 in 
esu (these values may be typical for the pulsar polar cap) , 
I i?i I is needed to be larger than 5 x 10 3 in esu to change 
the distribution of pairs. Even if \E\\ < \Eo\, the excited 
disturbances may change the 0-th distribution of pairs. 
However, we need numerical simulations, which can deal 
with non-linear process, in order to check how the reaction 
force modifies the distribution function fo- 

As a first step to deal with behaviour of pairs in an 
electric field, in this paper we have assumed that the back- 
ground charge distribution cancels out the modification of 
Eo due to injected electron-positron pairs. Of course, this 
simplification may not be appropriate for pulsars, while 
it makes computation easier. Inhomogeneous electric field 
might play an important role in plasma instability. Let 
us check whether our model can be used when the back- 
ground charge density is constant for s > sq. The charge 
density changes for s > so owing to the pair injection 
and electric field should be modified. The charge density 
decreases with distance as oc^ nos for s > So. In this 
case the variation of AEq over the moving distance of in- 
jected positrons before they turn back I^uq is ~ (Ruq) 2 Eq. 
Therefore, in the case of stable plasma (Ruo <C 1) the con- 
stant electric field is a good approximation even for the 
constant background. Although our simulations show a 
new possibility of plasma instability around pulsars, we 
need to simulate with an inhomogeneous electric field for 
Ruq ~ 1 in order to conclude whether instability occurs 
in actual situations on pulsars. In any case, the condition 
Ruo >~ 0.1 is an necessary (but not sufficient so far) 
condition to induce the instability. 

Let us consider implications for the pulsar polar cap. 
We suppose that the primary electron-beam is accelerated 
from z — 0, and its Lorentz factor becomes T at the pair 
production front (PPF) (z = L). In this case the aver- 
age electric field is \Eq\ — Tm e c 2 / eL. The braking time is 
expressed as ts = m c c/\eEo\ = L/cT. By curvature radi- 
ation an electron of T emits 2e 2 T/3rTi = M photons per 
unit time, where r is the radius of curvature of the field 
line. We express the number density of the primary beam 
as rib = fbriGJ- Assuming curvature gamma-rays imme- 
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diately turn into pairs, the pair-injection rate is approx- 
imated as ho ~ MfbTiQj. Adopting the average electric 
field, the ratio R becomes 



R = Wp^E — 



Uf h e 3 nBL 3 
3m c 4 frrr 2 



9 y in-5 f l/2 T -l/2 R l/2.3/2 -l/2 r -l 



(26) 
(27) 



where Tq = T/10 6 , and To. 3, Byi, L4, and r-j are in units 
of 0.3 s, 10 12 G, 10 4 cm, and 10 7 cm, respectively. The 
average Lorentz factor of pairs is at most KT 3 / (2m c cr) ~ 
2T|rf 1 . Therefore, we obtain 



Ru ~ 4 x 10- 5 fl /2 T-j /2 Bll 2 Ll /2 r- 3/2 Tl. 



(28) 



Even for the large values of T and L we can suppose in pul- 
sar models, we may expect 10~ 4 < Ruq < 10~ 2 at most. 
The electric field in the polar cap may be too strong to 
induce the instability, compared to the pair-plasma den- 
sity, while the approximation of constant Eq may be not 
so wrong in such cases. If the electric field is much smaller 
than the averag e one Trrip.c 2 leL at PPF as in the model 
in Aron' s group JFawlev et al. Ill977t Scharlemann et al. L 



Il978t lArons and Scharlemann L Il97£ ), Ruo can be large 
enough to induce the instability. In such models, however, 
the pair polarization is not so important to achieve the 
screening, while there is possibility that the instability af- 
fect the radio emission process. 

At present it is not clear that the electrostatic insta- 
bility we have considered in this paper is an important 
process for the screening of electric field above the pulsar 
polar cap. However, there may be extreme environments 
(magnetars etc.), where Ruq is as large as one. We expect 
that studies on the plasma instability in electric fields lead 
to opening a new approach to high-energy astrophysics. 
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Fig. 1. Distribution function of the toy model in §2. 
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Fig. 2. Distributions of electrons and positrons in the 
phase space. Electrons (positrons) are in the region u > uo 
(u < uo). The electric field £ (schematically shown as 
the central corrugation in the figure) generates the distur- 
bances F at the injection of pairs (u = u ), u = um, and 
u = u m . The arrows schematically show the transports of 
F (see text). 




Fig. 3. Electrostatic waves for RUN1 for r = 25, 165, and 
195. 




Fig. 4. Electrostatic waves for RUN1 around r = 200. 
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Fig. 5. Time evolution of Em for RUN1. 




Fig. 6. Charge-density distributions for r = 195 in RUN1. 
The dashed line is the electric field at that time. 




Fig. 7. Electrostatic waves for RUN2. The thick lines are 
for r = 80 and 90. The thin lines are for r = 330 and 340. 
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Fig. 8. Backward waves for RUN3. 




Fig. 10. Charge-density distributions for r = 300 in 
RUN3. The dashed line is the electric field at that time. 




Fig. 11. The 'average force' G(u) for r = 195 in RUN1. 
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R 


wo 


SM 


Si 


fci 


kA 


RUN1 


0.1 


10 


300 


120 


0.1 


1.0 


RUN2 


0.1 


300 


1500 


500 


0.1 


0.2 


RUN3 


0.1 


3 


500 


420 


1.0 


0.2 


RUN4 


0.01 


100 


1500 


200 


0.01 


1 


RUN5 


0.01 


30 


1500 


1320 


0.05 


0.01 


RUN6 


0.01 


3 


2000 


1820 


0.05 


0.01 



Table 1. Parameters and initial conditions in computa- 
tion. 





Ruo 


Type 


HUJp 


A/Z p 


Vph 


RUN1 


1.0 


A 


5 


6 


-2.8 


RUN3 


0.3 


C 


10 


0.9 


-1.1 


RUN4 


1.0 


A 


5 


8-10 


-2.8 


RUN5 


0.3 


C 


10-30 


1 


-1.1 



Table 2. Rough values of the growing time, wavelength, 
and phase velocity. The characters "A" and "C" represent 
the absolute instability and convective instability, respec- 
tively. 



